Granular Rayleigh-Taylor Instability: Experiments and Simulations 
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A granular instability driven by gravity is studied experimentally and numerically The instability arises as 
grains fall in a closed Hele-Shaw cell where a layer of dense granular material is positioned above a layer of 
air. The initially flat front defined by the grains subsequently develops into a pattern of falling granular fingers 
separated by rising bubbles of air. A transient coarsening of the front is observed right from the start by a finger 
merging process. The coarsening is later stabilized by new fingers growing from the center of the rising bubbles. 
The structures are quantified by means of Fourier analysis and quantitative agreement between experiment and 
computation is shown. This analysis also reveals scale invariance of the flow structures under overall change of 
spatial scale. 

PACS numbers: 45.70.Qj, 47.20.Ma, 89.75.Da, 47.11.-j 



Improved understanding of granular flows would be of es- 
sential importance to a range of industrial applications, to the 
study of geological pattern forming processes, and, in general, 
to the theoretical description of disordered media. 

As grains become smaller the effect of the interstitial fluid 
becomes more important. The result is a combination of dry 
granular dynamics and the hydrodynamics of the fluid. These 
systems give rise to a variety of exotic and most often poorly 
understood phenomena such as fluidization [1] and bubble in- 
stabilities [2], quicksand and jet formation [3], and sandwich 
structures in systems where different particle types segregate 
[4]. While the study of dry granular media has been exten- 
sive over the past decades, the exploration of fluid-granular 
systems has been of more limited scope. 

In the present Letter we study a granular analog of the 
Rayleigh-Taylor instability [5] in the sense that an interface 
instability arises as a heavier phase (the grains) displaces a 
lighter phase (the air). The experimental setup consists of a 
closed Hele-Shaw cell that confines air and fine grains. When 
the cell is turned upside down we observe the evolution of an 
initially sharp front formed by the falling grains. This evo- 
lution has three stages: (1) An initial decompaction phase is 
followed by (2) the formation of vertical falling fingers (the 
dark filaments in Fig. 1) organizing into cusp-shaped struc- 
tures that subsequently develop into (3) coarser finger-bubble 
structures. The last structures, seen in Fig. 1, represent a qua- 
sisteady state where two competing mechanisms produce a 
characteristic wavelength. The mechanism producing coarser 
scales originates as smaller bubbles lag behind bigger bubbles, 
giving rise to a finger merging process shaped like an inverted 
Y (see Fig. 1). This process resembles the coarsening seen in 
crystal growth [6] . The other mechanism, that produces finer 
scales and is active right from the start, is reminiscent of the 
tip splitting process seen in viscous fingering. It is manifested 
as thin filaments forming in the centre of the rising bubbles. 

Over the past few years a wide range of granular instabili- 
ties where various structures form along fluid-grain interfaces 
have been reported [2, 7]. Notably, the patterns formed by 
grains falling in a highly viscous liquid were investigated ex- 



perimentally and theoretically by Voltz et al. [8]. However, 
while the instability reported by Voltz shares its main quali- 
tative characteristics with the classical Rayleigh-Taylor insta- 
bility, i.e., a single dominating wavelength growing right from 
the start, our gas-grain instability grows through coarsening 
cusp- structures. 

The evolving structures further exhibit scale invariance un- 
der change of particle size, a feature which is supported both 
by observations and theoretical considerations. The simula- 
tions and experiments that are employed to shed light on the 
phenomenon at hand agree qualitatively, and to a significant 
extent, quantitatively, even though the model neglects both 
granular friction and the spatial direction normal to the Hele- 
Shaw cell. 

In the experiment a Hele-Shaw cell of inner dimensions 56 
mm x 86 mm x 1 mm is partially filled with air and monodis- 
perse beads (mass density 1.05 g/cm 3 , diameter 140 um) at at- 
mospheric pressure. The cell is rotated manually in about 0.2 
seconds from a lower to an upper vertical position to rapidly 
bring the layer of beads above the layer of air. Images of 
the evolving instability are recorded at a rate of 500 frames 
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Figure 1: (a) Experimental image and (b) numerical snapshot of a 
vertical Hele-Shaw cell where polystyrene beads (in black) of 140 
um in diameter displace air (in white). The cells are 56 mm wide and 
were rotated 0.2 seconds ago. 
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per second by a high speed digital camera with a resolution 
of 512x512 pixels; see Fig. 1(a). A simultaneous numerical 
snapshot is given in Fig. 1(b). 

The numerical model combines a continuum description of 
the air with a discrete description of the granular phase [9] . 
The effect of the granular phase on the air pressure is that of 
a deformable porous medium locally defined by the granular 
packing. The granular phase is modeled as discrete particles 
from which coarse grained solid fraction p(x,y) and velocity 
fields u(x,y) are obtained by means of a linear smoothing 
function [9] . This function distributes the mass and velocity 
of a particle among its four neighboring grid nodes (2.5 grain 
diameters apart). The continuum gas phase is described solely 
by its pressure P(x, y). The inertia of the gas, and hence its 
velocity field, is not considered. This is justified for small 
particle Reynolds numbers which is the case for our system. 

The pressure is governed by the equation [9] 

+ u • V?) = V • (P^VP) - PV • u , (1) 

where (j) = 1 — p is the porosity, n the permeability, u the 
granular velocity field, and p the gas viscosity. This equa- 
tion is derived from the continuity of air and grain mass, and 
Darcy's law with permeability k. The Carman-Kozeny rela- 
tion is assumed for the permeability, and the isothermal ideal 
gas law for the air. 

The grains are governed by Newton's second law 

dsr ^ VVP 

m— = rag + Fi , (2) 

at p 

where m, v, and V are respectively the mass, velocity, and 
volume of the grain. Contact dynamics [10] is used to cal- 
culate the interparticle force Fi which keeps the grains from 
overlapping. The dynamics of the grains are simplified by ne- 
glecting particle-particle and particle- wall friction. A lower 
cutoff is imposed on the solid fraction because the Carman- 
Kozeny relation is not valid as the solid fraction drops below 
0.25 [11]. This cutoff causes the permeability of the most di- 
lute regions of the system to be slightly lower than the true 
permeability. The effect is a slight overestimation of the pres- 
sure forces acting on the grains in the dilute regions. 

The spatiotemporal evolution of the air-grain interface in 
the experiment and simulation is presented in Figs. 2(b) and 
2(e), respectively. For every horizontal position x the interface 
height y(x) is defined in the following way: Moving down 
from the top, y(x) is the height where p(x, y) drops below a 
given threshold value [see Figs. 2(a) and 2(d)]. Because of the 
large density contrasts in the system the interface position is 
rather insensitive to the exact value of this threshold. For the 
experimental data the threshold is set on the gray levels in the 
images. 

The shape of the initial interfaces in Figs. 2(b) and 2(e) is 
quite different. The initial experimental interface has noise 
on all wavelengths, whereas the initial numerical interface is 
virtually flat with noise dominantly at smaller wavelengths. 



Perturbations introduced in the cell by the rotation and sud- 
den stop disturb the initial experimental interface. However, 
as the instability evolves the discrepancy between experiment 
and simulation reduces and the later interfaces are in better 
agreement. 

In order to give a more quantitative comparison of the inter- 
faces, the discrete Fourier transform with a Hamming window 
is applied on every second interface in Figs. 2(b) and 2(e) to 
produce the power spectra presented in Figs. 2(c) and 2(f). 
The power spectra are colored as their corresponding inter- 
faces, and the location of the maximum wave number for each 
power spectrum is indicated by a circle. While the maximum 
wave number of the numerical interfaces moves from high 
values to low values, the maximum wave number of the ex- 
perimental interfaces hardly moves at all, most likely because 
the experiment does not evolve from an initially flat interface. 
However, the experimental and numerical power spectra con- 
verge to approximately the same form when normalized. 

To study the coarsening of the observed structures quanti- 
tatively we perform an average over the solid fraction for the 
entire system: The discrete Fourier transform and the power 
spectrum of each horizontal line of p(x, y) is calculated. The 
averaged power spectrum, S(k), is then obtained by averag- 
ing over all these horizontal power spectra. An average wave 
number is defined as (k) = ^2 k S(k) • k/^2 k S(k) 9 where 
1/k is the wavelength. Likewise, we define the squared stan- 
dard deviation a 2 = ^ k S(k) -k 2 /^ k S(k) - (k) 2 . For the 
experimental data the image pixel values are used to estimate 
the solid fraction. 

Figure 3 shows the temporal evolution of (k) and a k (inset) 
for the numerical and experimental data. An additional set of 
experimental data is added to the plot. The numerical curve 
starts out with a significantly higher wave number than the 
experimental curves. However, the numerical data decreases 
monotonously until it coincides with the experimental data at 
about 0.2 seconds, after which the simulation and experiments 
show a similar coarsening behavior. Fingers are not observed 
in the experiment until 0.06 seconds have elapsed. During 
this time the grains merely form a dilute sheet that appears 
homogeneous on the experimental images. This particular ex- 
perimental initial state is caused by the sudden stop of the cell 
and is the most probable reason for the initial discrepancy be- 
tween simulation and experiment in Fig. 3. The fluctuations 
of (k) and a k are associated with the continuous nucleation 
and merging of fingers. 

We further investigate the behavior of the system as the 
overall spatial scale is changed: Keeping all length ratios and 
the particle number fixed, the size of the system will scale ac- 
cording to the diameter d of the grains. We measure the char- 
acteristic inverse length scale (k) as d is changed and observe 
a scale invariance of the evolution. A series of seven simula- 
tions are performed where d varies from 70 um to 490 um in 
steps of 70 um. The dimension of the numerical cell confining 
grains of 70 um in diameter is 28 mm x 34 mm. In these sim- 
ulations we have introduced the larger density of glass, rather 
than polystyrene, in order to minimize the numerical artifacts 
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Figure 2: (Color online) (a),(d) Interfaces superposed on the snapshots (t = 0.20 s) from Fig. 1. (b),(e) Spatiotemporal evolution of the air-grain 
interface from t = 0.00 s (bottom) to t = 0.29 s (top). (c),(f) Temporal evolution of the power spectrum of the interfaces. The circles indicate 
the location of the maximum wave number. 



associated with the solid fraction cutoff in the permeability. 
To compare, a series of experiments using polystyrene beads 
of 80, 140, 230, and 570 um in diameter, confined in Hele- 
Shaw cells that scale proportionally with d in all directions, 
are performed. 

Data-collapse plots of the rescaled mean wave number d(k) 
are shown in Figs. 4(a) (simulation) and 4(b) (experiment). 
These plots indicate that the characteristic size of the struc- 
tures is invariant when size is measured in units of d; the num- 
ber of grains that spans the width of the bubbles is the same 
for a wide range of grain sizes. 

Theoretically, the scale invariance of the product d(k) may 
be interpreted as follows: Compared to the other terms of 
Eqs. (1) and (2) the mdv/dt, F\ and PV • u terms may be 
shown to be small [12]. For that reason, these equations ex- 
hibit an approximate invariance under system size scaling. If 
we take SP to be the pressure deviation from the background 
pressure, express the velocity of grain i as v^ = 5vi + uo 
and the locally averaged granular velocity as u = Su + uo, 
where uo is the constant sedimentation velocity of a close 
packed system, this scaling may be expressed as x — > Ax, 
5P — > XSP, uo — > A 2 uo, Su — > XSu and k — > A 2 ft, where 
A is a scale factor. The structure formation of the system is 
governed by 5u and, since this velocity scales the same way 



with A as the length scales themselves, the evolution of any 
structure measured in units of d will be scale invariant. In par- 
ticular this is true for the structures measured by the length 
l/(k), and so d(k) is scale invariant. However, the invariance 
deteriorates both when particle size is increased, and when it 
is decreased. In the first case, the relative effect of granular 
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Figure 3: (Color online) Mean wave number (k) and standard de- 
viation (jfc (inset) for two experiments and one simulation, all using 
polystyrene beads of 140 um in diameter. 
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Figure 4: (Color online) Data-collapse plot of d(k) for a series of (a) 
simulations using glass beads, and (b) experiments using polystyrene 
beads. The grain diameters d are given in the legend box. The inset 
shows the evolution of (k). 



inertia is increased, in the second, the relative effect of the 
P V • u term is increased. 

The convergence of the numerical data-collapse in Fig. 4(a) 
is quite good. The deviation of the 70 um curve for small t is 
probably explained by the increase in the relative importance 
of the PV • u-term. The divergences of the 350, 420, and 
490 um curves for greater t in the same plot arise because the 
bubbles in the coarser packings disappear before they reach 
the surface due to the increase of u with A 2 [12]. The ex- 
perimental data in Fig. 4(b) have a wider distribution but col- 
lapses satisfactorily given the standard deviation error bars. 
The experimental data are obtained by averaging over three 
experiments for each diameter d. The standard deviation is 
calculated over a time window of 0.3 seconds. The accuracy 
of the experiments is at its lowest during the initial coarsen- 



ing of the structures, but as the mean wave number stabilizes 
around 0.2 seconds the accuracy improves. Nevertheless, the 
data points are, with a few exceptions, within a distance of one 
standard deviation from one another. The loss of precision for 
small times is most likely caused by the inaccuracy involved 
with the manual rotation. 



In conclusion, we have presented experimental and numer- 
ical results of a gravity-driven granular flow instability which 
is significantly different from its classical hydrodynamic ana- 
log. The simulations reproduce the characteristic shape and 
size of the experimentally observed structures and provide fine 
patterns in the early phase of the process that are not resolved 
experimentally. Data-collapse plots of the mean wave num- 
ber (k) indicate that the flow and the resulting structures are 
invariant when measured on a scale proportional to the grain 
diameter d for a range of diameters that spans from 70 um to 
570 um. 



* Electronic address: janlv@fys.uio.no 

[1] H. J. Herrmann, J.-R Hovi, and S. Luding, eds., Physics of 
Dry Granular Media (Kluwer Academic Publishers, Dordrecht, 
1998); D. Gidaspow, Multiphase Flow and Fluidization (Aca- 
demic Press, San Diego, 1994). 

[2] D. Gendron, H. Troadec, K. J. Mal0y, and E. G. Flekk0y, Phys. 
Rev. E 64, 021509 (2001); E. G. Flekk0y, S. McNamara, K. J. 
Mal0y, and D. Gendron, Phys. Rev. Lett. 87, 134302 (2001). 

[3] D. Lohse, R. Rauhe, R. Bergmann, and D. van der Meer, Nature 
432, 689 (2004). 

[4] C. Zeilstra, M. A. van der Hoef, and J. A. M. Kuipers, Phys. 
Rev. E 74, 010302 (2006). 

[5] Lord Rayleigh, Proc. London Math. Soc. 14, 170 (1883); 
G. Taylor, Proc. Roy. Soc. A 201, 192 (1950). 

[6] J. S. Langer, Rev. Mod. Phys. 52, 1 (1980). 

[7] 0. Johnsen, R. Toussaint, K. J. Mal0y, and E. G. Flekk0y, Phys. 
Rev. E 74, 011301 (2006); C. Chevalier, M. B. Amar, D. Bonn, 
and A. Lindner, J. Fluid Mech. 552, 83 (2006); F. Malloggi, 
J. Lanuza, B. Andreotti, and E. Clement, Europhys. Lett. 75, 
825 (2006). 

[8] C. Voltz, W. Pesch, and I. Rehberg, Phys. Rev. E 65, 011404 
(2001). 

[9] S. McNamara, E. G. Flekk0y, and K. J. Mal0y, Phys. Rev. E 61, 
4054 (2000); D.-V. Anghel, M. Strauss, S. McNamara, E. G. 
Flekk0y, and K. J. Mal0y, Phys. Rev. E 74, 029906(E) (2006). 
[10] F. Radjai, M. Jean, J. -J. Moreau, and S. Roux, Phys. Rev. Lett. 
77, 274 (1996). 

[11] A. A. Zick and G. M. Homsy, J. Fluid Mech. 115, 13 (1982). 
[12] J. L. Vinningland, 0. Johnsen, E. G. Flekk0y, R. Toussaint, and 
K. J. Mal0y (unpublished). 



